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ABSTRACT 

In this paper, we introduce the notion of subcell resolution, which is 
based on the observation that unlike point values, cell-averages of a dis- 
continuous piecewise-smooth function contain information about the exact loca- 
tion of the discontinuity within the cell. Using this observation we design 
an essentially non-oscillatory (ENO) reconstruction technique which is exact 
for cell averages of discontinuous piecewise-polynomial functions of the 
appropriate degree. Later on we incorporate this new reconstruction technique 
into Godunov-type schemes in order to produce a modification of the ENO 
schemes which prevents the smearing of contact discontinuities. 
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1. INTRODUCTION 


In [7], [8], and [9] we have introduced a class of essentially non- 
oscillatory (ENO) schemes that generalizes Godunov's scheme [2] and its second 
order extensions ([10], [1]) to high order of accuracy. 

In this paper we present a modification of the ENO schemes which is 
designed to prevent smearing of linear discontinuities. 

Let { lj x tt n ,t n+1 )}, where Ij » [ x j-l/2» x j+l/2l » x a = ah > \ = kx » be 
a partition of R xR + . Let u? be the the "cell-average" of u at time 
t n > i* e. , 


( 1 . 1 ) 


— n 1 


u j ’ h ( 

J 


The cell-average of the solution to the initial value problem 


( 1 . 2 ) 


u t + f(u) x = 0, u(x,0) = uq(x) 


satisfies 


(1.3a) 


“f 1 " - 1[?(x J +l/2' t n iu) ‘ ?(x j-l/2- t n ;u)1 


where X = t/h and 


(1.3b) 


1 T 

f(x,t;u) = — / f(u(x,t+n))dn> 

T r\ 


The ENO schemes can be written in the standard conservation form 


- 2 - 


(1.4a) 


n+1 n /-£ 
v. = v ; - X(f, 


r^ (f j+ i/2- f ri/2 )E tE h (T) ‘ v V 


here denotes the numerical solution operator and ^j+1/2* tbe 


irical flux, denotes a function of 2k variables 


(1.4b) 


f j+l/2 f(v j-k+l’ ***’ V j+k^ 


which is consistent with the flux f(u) in (1.2), in the sense that 
f(u,u,»««,u) = f(u). Unlike standard difference schemes, in the ENO 

schemes is a high-order approximation to the cell-average u?, and not to 
the point value u(xj,t n ). Setting v” = u^ in the numerical scheme (1.4) 
and comparing it to relation (1.3), we see that if the numerical flux 
fj + l /2 * ,u^ +k > can be expanded as 

(1.5a) T(u^_ k+1 ,... ,u^ +k ) =1/ f(u(x j+1/2 ,t n +n))dn + d(x J+1/2 )h r + 0(h r+1 ) 


then the truncation error 


(1.5b) u? +1 - [E k » u^Jj = \[d(Xj +1 ^ 2 ) - d(Xj_^ 2 )]h r + 0(h r+1 ), 

is 0(h r+ *) wherever d(x) is Lipschitz continuous, i.e., the scheme (1.4) 
is r~th order accurate in the sense of cell averages. 

The most important ingredient in the ENO schemes is a procedure to re- 
construct a piecewise-smooth function w(x) from its given cell-averages 
(Wj) . This reconstruction, which we denote by R(x;w), is a piecewise- 
polynomial function of x that has a uniform polynomial degree (r - 1) and 


— ■*— 


satisfies: 

(i) At all points x for which there is a neighborhood where w is 
smooth 

(1.6a) R(x;w) = w(x) + e(x)h r + CKh 1 ^*); 

(ii) Conservation in the sense of 

1 V 1 / 2 

(1.6b) ^ / R(Xj + C;w)d£ = w^; 

X j —1/2 

(iii) It is essentially non-oscillatory 

1 4-P 

(1.6c) TV(R(.;w)) < TV (w) + 0(h ), P > 0 


where TV denotes total variation in x. 

Using the reconstruction (1.6) we can express the abstract form of the 
ENO schemes by 


(1.7a) 


TT, (t)»w = A(I)»E(r )»R(» ;w). 
n 


Here A(I) is the cell-averaging operator 


(1.7b) 


A(I).w 


VT 


f w(y)dy 
I 


and E(t) is the exact evolution operator of the IVP (1.2), i.e., 
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(1.7c) u( • ,t) = E(t)u 0 . 

We note that (1.7a) with the piecewise constant reconstruction 
(1.8) R(x;w) = Wj for Xj_^ 2 <. x <. x j+i /2 

is exactly the first order accurate Godunov's scheme [1]; (1.7a) with the 

piecewise linear reconstruction 

(1.9a) R(x;w) = w ^ + S^(x - x ^ ) for x j_i /2 1 x < x j+i/2 

where 


(1.9b) S. = w (x.) + 0(h), 

J x J 

is the abstract form of the second order accurate extensions to Godunov's 
scheme described in [10], [1], and [7], 

In the first order case (1.8), the scheme (1.7a) can be expressed in the 
conservation form (1.4) with the numerical flux 


( 1 . 10 ) 


- 2 - r R, n n » 

Vi/2 • f <v j-Vl ); 


here f R (u^,U 2 ) is an approximation to the flux at the origin in a Riemann 
problem with to the left and U 2 to the right. 

In the second order case (1,9), the numerical flux of the abstract scheme 
cannot be expressed in a simple closed form, and we approximate it by 
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(1.11a) 


T = ■c® N0 = f R, L R . 

j+1/2 j+1/2 ^ V j+l/2’ v j+l/2^ 


(1.11b) 


j+1/2 


= v" + h(l - 


Xa? )Sj/2, 


V j+l/2 v j+l ” h( - 1 + Xa j+l^ S j+l^ 2; 


here a a = f'(v a ). 


In this paper we pay special attention to the second order accurate 
scheme (1.11), because at present this seems to be the state of the art. This 
class of second order schemes (with various choices of S a ) performs rather 
well in smooth regions and shocks. However, it exhibits excessive smearing of 
linear discontinuities, i.e., contact discontinuities. Usually such discon- 
tinuities are smeared more and more in time at the rate 0(n*^), where n 
is the number of time-steps. To understand this smearing we note in (1.7) 
that whenever a discontinuity in the reconstruction R is propagated by the 
evolution operator E into the interior of the cell, then the cell-averaging 
operator A(I) replaces this sharp discontinuity by a smeared transition. In 
the linear case there is nothing to stop this process and therefore it goes on 
forever. In the case of a shock wave, the fact that the characteristics con- 
verge into the shock counteracts the smearing, and a steady progressing 
profile is obtained. 

The above observation is the basis for the artificial compression concept 
[3], In order to prevent the excessive smearing of a linear discontinuity one 
can artificially induce convergence of the numerical characteristic field at 
each monotone strip of the solution. This can be accomplished by modifying 
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the expression of the slopes S^, or by adding a corrective term to the 
numerical flux (1.11) (see [4]). The main advantage of artificial compression 
is that it is easy to use. The primary disadvantage is that one has to be 
extra careful (which also means to do a lot of checking...) not to generate 
unphysical discontinuities by applying it too strongly, where it need not be 
applied at all. We refer the reader to [6] for more details. 

The piecewise-parabolic method (PPM) of Colella and Woodward [1] includes 
a mechanism to detect contact discontinuities and to correct the scheme by 
using a "steeper" reconstruction. The PPM proved itself to be a robust high 
resolution scheme in a large number of numerical tests [12]. In this paper we 
present a technique, which we call "subcell resolution," that is close in 
spirit to the PPM but is somewhat different in its methodology. 

The present scheme is a "souped up" version of (1.11) in which the linear 
advection part is boosted to third-order accuracy (in Lj-sense) and is 
capable of propagating linear discontinuities perfectly (within 3rd order 
accuracy). The main ingredient in the new method is the observation that the 
information in cell-averages of a discontinuous function, unlike that of point 
values, contains the location of the discontinuity within the cell, e.g. , the 

(“L 3 i • 1 
"j -{“m 3 ■ 0 


cell-averages w^ 


(1.12a) 


with u M between 


u^ and u^, are identical to those of the step-function 
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Using this observation we can modify the ENO reconstruction of [8] to recover 
exactly any discontinuous quadratic function from its cell-averages. 

In order to retain the relative simplicity of the numerical scheme (1.11) 
we use the new reconstruction to correct only the linear advection part. The 
new numerical flux is 


(1.13) 


T _ -rENO - 

j+1/2 " j+1/2 g j+l/2 ; 


here §j+i/2 t * ie ^ ux through Xj+ 1/2 due to the linear advection of 
the difference between the modified reconstruction and the piecewise-linear 
one (1.9). In the constant coefficient case the scheme (1.13) is exact for 
discontinuous quadratic initial data. 

Later on in this paper we present the extension of the "subcell resolu- 
tion" concept to any finite order of accuracy, and also extend the scheme to 
the Euler equations of gas dynamics. 


2. ENO RECONSTRUCTION 

In this section, we describe one of the techniques to obtain an ENO re- 
construction. Given cell-averages {w.} of a piecewise smooth function 
w(x), we observe that 
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__ x j +1/2 

(2.1a) h W j = f w(y)dy = w ( x j +1 / 2 ) - VKx^j^) 

X j — 1/2 

where 


x 

(2.1b) W(x) = / w(y)dy 

X 0 

is the primitive function of w(x). Hence we can easily compute the point 
values t w ^ x i+i/2^ by summation 


(2.1c) 


W(x i+l/2 ) 


i 

h I 

j=i 


0 


w . . 

3 


Let H m (x;u) 
accurate to order 


be an interpolation of 
rn j i • e • | 


u 


at the points 


{yj}, which is 


(2.2a) H (y ;u) = u(y ), 

m 3 J 

(2.3b) H (x;u) = ii- u(x) + 0(h" H ' 1 " £ ), 0 < SL < m. 

, Jc m . £ — — 

dx dx 

We obtain our "reconstruction via primitive function" technique by 
defining 

(2.4) R(x;w) = H^(x;W). 

Relation (1.6a) follows immediately from (2.3b) with l ~ 1 and the defini- 
tion (2.1), i.e.. 



Relation (1.6b) is a direct consequence of (2.3a) and (2.2), i.e., 


1 X 3 +1/2 d 

A(I )R(. ;w) = H r (x,W)dx 

X j — 1/2 


= ¥ [H r (x j+l/2 jW) " V X j-l/ 2 ;W] = ¥ tW(x j+l/2 ) “ W(x j-l/2 )] = W j* 

To obtain an ENO reconstruction, we take H r in (2.4) to be the new ENO 
interpolation technique of the author [5], In this case, H m (x;u) is a 
piecewise-polynomial function of x of degree m, which is defined (omitting 
the u dependence) by 


(2.5a) 


H m ( X ;u) = q j+1/2 (x) for y . < y < y j+1 


where ^j+ 1/2 * s t ^ ie un ^- < l ue polynomial of degree m that interpolates u 
at the m+1 points 


(2.5b) 


S m (1) 5 


for a particular choice of i = i(j) (to be described in the following). To 
satisfy (2.3a), we need 
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Vl/2 (y j > " “'V’ q j+l/2 ( Vl > ‘ u(y j+l ,; 

therefore, we limit our choice of i(j) to 

(2.5c) j - m+1 < i(j) < j. 

The ENO interpolation technique is nonlinear: At each interval 

[yj,yj +1 ], we consider the m possible choices of stencils (2.5b) subject to 
the restriction (2.5c), and assign to this interval the stencil in which u 
is "smoothest" in some sense; this is done by specifying i(j) in (2.5b). 

The information about the smoothness of u can be extracted from a table 
of divided differences. The k-th divided difference of u 

(2.6a) u[y i» y i+r*’* ^i+k 1 H u[S k (i)] 

is defined inductively by 

(2.6b) u[S Q (i) ] = u(y t ) 

and 

(2.6c) u[S k (i)] = (u[S k _ 1 (i+1) ] - u[S k _ 1 (i)])/(y 1+k -y i ). 

If u(x) is m times differentiable in [y^y^.^] then 
(2.7a) u[S m (i)] =J-j-i/ m \s), for some 5 < y i+ni * 
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If u^P^Cx) has a jump discontinuity in [yi>yi+nJ then 

(2.7b) u[S m (i)] = 0(h~ m+ ^ ) [u^^]), 0 <p < m-1 

([u<P>] in the RHS of (2.7b) denotes the jump in the p-th derivative). 

Relations (2.7) show that |u[S m (i)]| is a measure of the smoothness 
of u in S m (i), and therefore can serve as a tool to compare the relative 
smoothness of u in various stencils. The simplest algorithm to assign 
S m (i(j)) to the interval [yj»yj+j] is the following: 

Algorithm I. Choose i(j) so that 

(2.8) | u [ S ( i ( j ) ) ] | = min { |u[S (i) ] |> . 

m j-m+l<i<j “ 

Clearly (2.8) selects the "smoothest" stencil, provided that h is 
sufficiently small. 

In order to make a sensible selection of stencil also in the "pre- 
asymptotic" case, we prefer to use the following hierarchial algorithm: 

Algorithm II. Let i^(j) be such that S| c (ij c ( j ) ) is our choice of a (k+1)- 
point stencil for [yj,yj+il» Obviously we have to set 


(2.9a) 


il(j) - j 
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To choose ik+l^^» we consider as candidates the two stencils 

(2 - 9b > ti ' WVJ> - »• 

<2.9c) S^ +1 - S k+1 (l k(J )), 

which are obtained by adding a point to the left of (or to the right of) 
Sj^ikXj)), respectively. We select the one in which u is relatively 

smoother, i.e., 

(l k (j)-l If |u[S b + ,J| < |u[S« +1 ]| 

(2. 9d) i k+ ,(j) = < 

(i^Cj) otherwise 

Finally we set i(j) = i m (j). 

Using Newton's form of interpolation, we see that the polynomials 
(q k (x)}, 1 £ k £ m, corresponding to the stencils S k = S k (i k (j)) selected 
by Algorithm II, satisfy the relation 

(2. 9e) q k+l^ = q k^ + u f sk+1 ^ 11 k ( x_ y)* 

This shows that the choice made in (2.9d) selects 9k+l to t ' ie one tkat 

deviates the least from q k . It is this property that makes Algorithm II 

meaningful also for h in the pre-asymptotic range. 
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3. ENO RECONSTRUCTION WITH SUBCELL RESOLUTION 

In this section, we show how to modify the ENO reconstruction of the pre- 
vious section so as to allow for the recovery of discontinuities in the in- 
terior of the cells. To illustrate the procedure, we first consider a dis- 
continuous piecewise polynomial function w(x) of the form 




P L (x) X < Xj 

(3.1a) 

w(x) = 1 

> 


1 

P R (x) X > x d 

where P^(x) 

and P R (x) are polynomials of degree less or equal s 

(3.1b) 

deg(P L ) £ s, deg(P R ) _< s. 


We assume that w(x) is actually discontinuous at x^, i.e., 

(3.1c) P L^ X d^ * P R^ X d^ 

and that the discontinuity is located in the interior of the interval Iq 
( 3. id) x_ 1/2 < x d < x 1/2 . 

(See Figure la.) 

Next we denote the cell-averages of w(x) in (3.1) by {w^} an£ * con “ 
sider the ENO reconstruction R(x;w) applied to these data. To simplify 
our presentation let us denote the polynomial defining R(x;w) in the 
cell Ij by Rj(x;w). Clearly, provided h is sufficiently small, the 
stencils assigned to cells {I j }, j ^ 0, are selected from the smooth part of 



the function. Therefore, it follows from (2.3) - (2.5) that 


(3.2) 


Rj(x;w) = P L (x) + 0(h r ) for j £ -1 

R.(x;w) = P p (x) + 0(h r ) for j > 1. 

J K — 


R(x;w) in Iq does not introduce spurious oscillations, however it does not 
provide an accurate approximation to the discontinous function w(x) either 
(see Fig. lb). Using (3.2) and the information contained in the cell- 
average Wq, we can easily rectify this situation as follows: We extend 

R_j(x;w) to a point z in Iq from the left, and extend R^(x;w) to 

z from the right and then approximate the location of the discontinuity in 
the cell Iq by finding a value of z that will fit the cell average 
Wq (see Fig. lc). This is done by finding a root of the algebraic equation 
Fq(z) = 0 where 

1 2 _ Xl / 2 

(3.3a) F Q (z) = {/ R -;1 (x;w)dx + / R J (x;w)dx} - w Q . 

X — 1/2 z 

When h is sufficiently small, the data near the cell Iq approach those of 
a step-function. Therefore, as in (1.12) we expect to have 


(3.3b) 




F 0 ( x 1/2 


) < 0 


and a single root of f q(z) =0 in Iq; we denote this root by 0 q 
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It follows from (3.2) that 


(3.3d) 


l e o “ x dl = 0(hr> * 


What we mean by "ENO reconstruction with subcell resolution" is the modified 
ENO reconstruction R(x;w) which is defined in this case by 


(3.4a) 

(3.4b) 


R.(x;w) = R. 

3 3 

(x;w) 

for j £ 0 

* _ 1 

R _i (x;w) 

for 

x -l/2 < x < 9 0 

II 

/-\ 

1* 

X 

P*p 

R^ (x;w) 

for 

e 0 < x < x l/2 


Clearly it follows from (3.2) and (3.3d) that R(x;w) is an 0(h r ) 
approximation to w(x) in the Lj sense, i.e., for any a and b 


(3.5) 


j |R(x;w) - w(x)|dx = 0(h r ). 
a 


We observe that if the polynomial degree s in (3. Id) is less or equal 
(r - 1) then the primitive functions of and Pj^ are polynomials of 
degree less or equal r, and therefore H r (x;w) in (2.5) is exact except at 
Iq. Hence, 


(3.6a) 

(3.6b) 


Rj(x;w) « P L (x) for 3 <, “1 


R.(x;w) = P D (x) for j > 1, 
J K “ 


9 o ■ x d 


and consequently 


in (3.3) 


Thus we have shown 
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(3.6c) 


s < r -1 R(x;w) = w(x). 


We turn now to describe the algorithm defining R(x;w) for a general 

Tn _ 

piecewise-smooth function w(x). As in the previous example we take R(x;w) 

in to be R.(x;w), unless 1^ is suspected of having a discontinuity 

J J J 

of w(x) in its interior. In the latter case we check whether 


(3.7a) 




) < 0 , 


where 


(3.7b) 


F j (z) 


•in 


h R j(x;w)dx + J 

X j-l/2 z 


l j+l/2 


Rj +1 (x;w)dx} 


- w . 


If (3.7a) holds, then there is a root z = 0^, 


(3.7c) 


W ■°’ 


X j — 1/2 < < * 3 + 1/2 


in the cell I j , and we define R(x;w) in this cell to be 

S R. !(x;w) for < x < 9 . 

_ 

R^ +1 (x;w) for 0j < x < x j +1 / 2 


If (3.7a) does not hold (which means that either there is no root in I j , 

a 

or that there is an even number of roots in Ij), we ta ^ e R^(x;w) to be 
R^(xjw). 

Let Oj be some measure of the non-smoothness of the reconstruction 
R(x;w) in I j , e.g. , 
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H k _ 

(3.8a) a. = |^-r- R(x, ;w) | , 1 < k < r-1, 

2 dx K 2 " 

or a combination of such derivatives. Our algorithm identifies cells which 
are suspect of harboring a discontinuity of w(x), as those which attain a 
local maximum of "non-smoothness" of the reconstruction, i.e.. 


(3.8b) 




and 


We summarize Che algorithm defining 


a j - Vi* 

A 

R(x;w) by: 


If 


(3.9a) 

a , >a. ,, 0 . >0.,, and 
j j-l j - J+l 


then 

. _ ( R j-i <x;:r) 

for x j_i/2 ^ x ^ 

(3.9b) 

R (x;w) = j __ 

2 VR J+1 (x;w) 

5 


for Qj < x < x j+1/2 

otherwise 



(3.9c) 

Rj(x;w) ■ 

= Rj(x;w). 


In an Appendix, we present analysis which motivates the choice of condi- 
tion (3.8b). In the following we make several remarks and observations 

a 

about R(x;w), the "ENO reconstruction with subcell resolution." 
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(1) R(x;w) is indeed essentially non-oscillatory (ENO). This follows from 
the fact that local maxima are isolated, i.e., if (3.8b) holds for I j , 

A _____ 

it cannot hold for neither *j-l nor lj+i» Consequently, if R(x;w) 
is defined in Ij by the discontinuous (3.9b), then in I and 
Ij + ^ it is defined by (3.9c), i.e., 

a __ A _ 

(3.10) Rj.^xjw) = Rj_ x (x;w) , R j+1 (x;w) = R^ +1 (x;w). 

(2) If, as in the example (3.1), there is a discontinuity of w(x) in the 
interior of Iq, then 

(3.11a) a Q = 0(h" k ), a_ { = 0(1), Oj = 0(1) 

where k is the order of the derivative used in (3.8a). 

Therefore, provided h is sufficiently small, we get that 
(3.11b) a Q > o-l, Oq > Op 

and also as in (3.3b) 

(3.11c) F 0 (x -1/2 ) * F 0 (x 1/2 ) < °* 

This shows that R(x;w), as defined by the general algorithm (3.9), will 


recover any real discontinuity of w(x). 
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(3) We observe that condition (3.9a) may hold also in the smooth part of 
w(x) near a local maximum of |d*S?/dx^|. In this case the algorithm 
places a discontinuity at 0j in the interior of Ij. However, 

because of the smoothness of w(x) there and (1.6a), we have 


(3.12) 


R j±i (e j ; ") = w( V + 0(hr) * 


Consequently the jump is of the size of the reconstruction error 0(h r ). 
We recall that the original ENO reconstruction is discontinuous at 


(x 


j+1/2 


}. Therefore, the effect of the algorithm (3.9) is to replace the 


two discontinuities at Xj-i /2 and Xj+j/2 a s i n g^ e one at 0j. 


A 

(4) We observe that in order to evaluate R^(x;w) in a cell Ij which 

contains a discontinuity at 0,. we have to find out whether x > 0. 

3' 3 

or x < 0 . Assuming 0. to be the only root of F.(0.) = 0 in 
3 3 3 3 

Ij, as is the case for a real discontinuity, we can use the logic of the 

interval-halving technique to evaluate R^(x;w) without actually 

computing 0.. To do so we calculate F-(x) and compare its sign with 
J ^ 

that of Fj( X j_ 1/2 ) or Fj(x j+1/2 ): 

Rj_ 1 (x;w) if Fj(x)«Fj(x^_ 1 ^ 2 ) > 0 


(3.13) 


R(x;w) « 


Rj + 1 (x;w) otherwise 


4. A SECOND ORDER ACCURATE ENO SCHEME WITH SUBCELL RESOLUTION 

In the following we describe how to incorporate the reconstruction with 
subcell resolution into the ENO schemes, so as to improve their resolution of 
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linear discontinuities. In this section, we present the derivation of (1.13), 
which is an improved version of the second order accurate MUSCL-type scheme 
(1.11). In the next section, we shall generalize these ideas to any order of 
accuracy. 

We start with the piecewise-parabolic reconstruction R(x;w) which is 
defined by (2.4) with r = 3, i.e., 

(4.1a) R(x;w) = ^ H 3 (x;W) 

where W is the primitive function (2.1b) of w(x), and H^xjW) is a 

piecewise-cubic ENO interpolation. Let i = i(j) be the left endpoint of the 

stencil ( x i-i/ 2 » x i+l/ 2 » x i+ 3 / 2 » x i+ 5 / 2 ^ assigned to the interval 
( x j-l/2> x j+l/2^ (2.8) or (2.9) or some other ENO technique; 

j - 2 i <_ j. The parabola describing R(x;w) in ( x j-i/ 2 > x j+l/2^ * s 
given by 

d 3 H 3 (x;W) _ _ _ 2 

(4.1b) C = 3 ( w i+2 ” 2 w i+i + w i)/ h 

j <j x J 

d 2 H~(x;w) _ __ 

(4.1c) Sj —5 | x=x = (w i+1 - w t )/h + (j - i - l/2)hCj 

QX J 

2 

(4. Id) R(x;w) = (Wj - C^) + S^(x “ x j) + y C j^ x ” x j) 2, 

a 

Using the algorithm described in (3.9) we now define R(x;w), which 
modifies R(x;w) in (4.1) so that it includes a discontinuity in the 
interior of each cell 1^ which meets condition (3.9a). Ideally, we would 
like to use the scheme 
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(4.2) Vj +1 = A(Ij)E(t)R(» ; v n ), 

which is third order accurate in L^-sense (but only second-order accurate in 
the maximum norm)* However, the proper approximation of the nuraerial flux of 

(4.2) is much more complicated than (1*11), since it is one order more 
accurate in time, and on top of it one has to account for discontinuities 
crossing the boundaries of the cell during the time-step. Bearing in mind 
that the main fault in the MUSCL-type scheme (1.11) that we want to correct is 
the smearing of linear discontinuities, we settle for the simpler second order 
scheme (1.13), which will be identical to (4.2) only in the constant co- 
efficient case. 

Our basic scheme remains 

(4.3a) Vj n+1 = A(Ij )E(t )L(» ;v°) 

with the piecewise linear reconstruction 


(4.3b) L(x;w) = w^ + S^(x - x^) x e I ^ ; 

as in (1.11) we approximate its numerical flux by 


(4.4a) 


_ -zENO f R, L R > 

j+1/2 " j+1/2 ' v j+l/2» v j+l/2'» 


where 


V j+l/2 


= v? + h(l - Xa")s"/2, 


j ' 3 ' 


J+1/2 


- VI - h(1 


+ Aa" )S? 




j+1 


/ 2 . 


(4.4b) 
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We take Sj in (4.4b) to be (4.1c), and observe that by (2.3b) 
(4.4c) Sj = w x (xj) + 0(h 2 ), 


which is one order higher than (1.9b). Consequently, as the UNO scheme of 
Harten and Osher [7], this scheme is truly second order in all L p norms 
(unlike TVD schemes which are first-order in L and second order only in 

) . 

We introduce subcell resolution into the scheme by modifying its numer- 
ical flux to be (1.13), i.e., we consider the scheme 


(4.5a) 


n+1 n , t~F t \ 

V J ' v j - X(f j+ l/2 ' f j-l/2 ) 


(4.5b) 


T = ^ENO “ 

j + 1/2 j+1/2 g j+l/2* 


In the constant coefficient case, 


(4.6) 


u. + au = 0, a = 
t x ’ 


= const. 


we define gj+j/2 to 1,6 


( 4 * 7a ) g j+l/2 t l [R(x j+l/2 


^ n \ 
at; v ) 


I.(Xj + i /2 - at;v n )]dt. 


Since in the constant coefficient case 


(4.7b) 


■^ENO a r ' * n Nj . 

1 j+1/2 “?/ L(x j+l/2 " at! v )dC 
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we get that 

(^•7e) f j+1/2 = ~ l ^ X j+l/2 " at; v ^ dt * 

This shows that in the constant coefficient case (4.6), the scheme (4.5) is 
identical to (4.2). Consequently, it is exact for discontinuous parabolic 
data of the form (3.1). 

Next we derive an expression for gj+ 1/2 constant coefficient 

case (4.7a); this will later be generalized to the nonlinear case by "freez- 
ing" the charasteristic speed within the cell. 

In the constant coefficient case, (4.7a) can be rewritten as 

. x j+l/2 

(4.8) % j+1/2 “ ~ / [R(y;v°) - L(y;v n )]dy. 

X j+l/2-ax 

First let us assume a > 0. When 

(4.9a) R(x;v n ) = Rj(x;v n ) in Ij, 

then 

i V1/2 n , 

(4.9b) gj+l/2 = 7 / [Rj(y»v ) - L(y;v )]dy = jj (v - 1 )(2v-l )ti ; 

X j+l/2-at 

here v = Aa. 

A __ 

When there is a discontinuity in the interior of Ij, Rj(x;v n ) is 


given by (3.9b), i.e.. 
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(4.10a) 


R(x;v n ) = 


Cx; v n > x j_i /2 < x < Qj 


R. . (x;v n ) 0 . < x < x 


j+1 


j - A - j+1/2 


and we have to find out whether (x jy 2 ~ ax ^ ^ ar £ er on smaller than 

0^. We recall that for (4.10a) to hold the relations 


(4.10b) 


F j^ 9 j^ ~ °* F j^ x j-l/2^* F j ^ X j+l/2^ — ° 


have to be satisfied, where F^(z) is defined by (3.7b). Using the basic 
idea of the interval-halving method for calculating a root of algebraic equa- 
tions, we can find out in which of the two cases we are in without actually 
calculating 0_. (see Remark (4) at the end of the previous section). All 
we have to do is to compute F j^ x j+l/2~ at) and compare its sign to that 
of F j( x j-i/2^» i* 6 ** 


(4.10c) F j^ x j+l/2 “ aT /?) > 0 


j v j-1/2' 


"j+1/2 


- at < 0 


(4. lOd) F j (x j+l/2 " ax),F j (x j-l/2 ) i 0 ^ x j+l/2 " aT > 9 j‘ 


To express the integral in (4.8) let us introduce the notation 


y 2 

(4.11) b m (y ] ,,y 2 ) = / R ffl (x;v n )dx 


K%-I2r C m^ + T S m^ + iV J l y ;-x 


1 „ 2 . 1 „ 3 1 y 2 -x m 


m 


In case x j + i/2~ ax 2. we S et: from (4 • 8 ) and (4.10a) that 


(4.12a) g j+1/2 = j- (b j+1 U j+1/2 ' aT > x j+l/2 ) " ax [v j + 7 (h " ax)S j ] >* 
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When x j+i/2 “ aT ^ 0 j> we use the fact that 


X j+l/2 . 


X j-l/2 


R(x;v n )dx = v” 


to express the integral in (4.8) by 


j+1/2 „ 


/ R(x;v n )dx = hv a - j 

v — rjT J v 


X j+l/2~ aT 


x j+l/2” aT 


X j-l/2 


R(x;v n )dx = hv a - / 


X j+l/2“ ax 


X j-l/2 


(x; v )dx; 


Rearranging terms we get in this case 


(4.12b) g j+1/2 -{(h - ax ) (Vj ~ -3— Sj ) - b j _ 1 (x j _ 1/2 , x j+1/2 - ax)}, 


To summarize, the definition of Sj+i/2 ln the case a > 0 is! 


(4. 13a) + 


s j+1/2 TZ 


(v — 1) (2v — 1 )h^C . , 


unless the discontinuity condition (3.9a) holds for Ij; in the latter case we 


define 


(4. 13b) + 


( h ' aTK Vf S j>" b j- 1 < x j- 1/2 "J+l/2' aI ) lf F j (x j.l/2' aT),F j (:, j-l/2 ) > 0 

b j +1 ( x j +1 / 2 -aT » x j+i/2^“ aT t v j + y(h-ax)Sj] otherwise 


Tg j+l/2 = 


Similarly for a < 0 we get: 
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(4.13a)~ t gj_i /2 = JJ (v+l)(2v+l)h 2 Cj , 

unless the discontinuity condition (3.9a) holds for I j ; in the latter case we 
define 


(4.13b)" 

! -bj_ 1 (x^_ 1/2 ,Xj_ 1/2 -aT )-at [Vj-j(h+aT)]Sj if F j (x j-i/ 2 )* F j (x j-i/2 -aT ) > 0 

n ax 

b j +1 ( x j_ 1 /2" aT » x j+i/2 ) “ (h+ax)(v j“ T~ S j ) otherwise 


Note that the expressions (4.13) are formulated as the contribution of the 
cell Ij to the numerical flux g. Thus, if a > 0, the contribution of 

A 

from the Ij cell goes to &j+i/2* wb ile if a < 0 this contribution goes 
to 

In section 6, we extend the subcell resolution ideas to the Euler equa- 
tions of gas dynamics. Since shocks are highly resolved by the original ENO 
scheme, we apply subcell resolution only to the linearly degenerate 
characteristic field in order to improve the resolution of contact dis- 

continuities. In this case the characteristic speed, which is the velocity of 
the flow, is not a constant but a function of the solution itself. Neverthe- 
less, we use the same expressions as in (4.13), except that a in 1^ is 

A 

replaced by a j , and v = Aa^. We compute the corrective flux g in the 

A 

following way: First we preset g * 0, and then we sweep over the mesh and 

collect contributions to g from each cell: If aj > 0 we add the RHS of 

j A 

(4.13) to Sj+i/2’ a j < 0 we add the RHS of (4.13)" to Sj-i/2* 

Note that if a j+i < 0 and a^ > 0, then g ^ +1/2 ® ets contr: *- bu tions from 

both Ij and Vl* 
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5. EXTENSION TO HIGH ORDER OF ACCURACY 

In this section, we describe the extension of the ENO scheme with subcell 
resolution to arbitrarily high order of accuracy. As in the second order case 
(4.5), we introduce subcell resolution to the high order accurate ENO schemes 
via a corrective flux ^j+i/2’ we consider the modified scheme 


(5.1a) 


n+1 

v j 


v j “ X(f j+l/2 “ f j-1/2^ 


(5.1b) 


T -gENO ~ 

j+1/2 j+1/2 g j+l/2* 


First, we describe briefly the derivation of ^j+1/2* re ^ er 

reader to [8] for more details. Let L(x;w) be an r-th order accurate 
reconstruction of w(x), such that 


(5.2a) r- f L(x;w)dx = w . 

h I. 3 

3 

As before we denote the (r-l)-th degree polynomial describing L(x;w) in 
I. by L.(x;w), and express it in the following Taylor expansion 


(5.2b) 


r-1 T3, 

L (x;w) = l 7 (x 
3 k=0 




,1_ L(X 
dx J 


Let Uj(x,t) 


be the solution to the initial value problem 


(5.3a) 


u + f (u) = 0, u(x,0) = L (x;w). 

t X J 


Since the initial data in (5.3a) are analytic, we know from the Cauchy 
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Kowalewski theorem that Uj(x,t) exists uniquely and it is analytic for some 

time 0 < t < t . Therefore, its Taylor expansion around x = x. and t = 0 
— — c J 


(5.3b) 


V*’ 0 ‘ j n rr L 0 V-. 

J k=0 m=0 


^.m, N k-m 
t (x-x.) , 


(5.3b) 


— u.(x ,0), 


m,k-m 3 t m 3 X k-ra j J 


is valid for 0 < t < t and x sufficiently close to Using (5.3) we 

— — c J 

define v.(x,t) to be the truncated Taylor expansion 


(5.4a) 


r-1 

v.(x,t) = l 
J k=0 


1 K 

7T l 

m=0 


( k ) 


m,k-m 


t m <x - 


V 


k-m 


(5.4b) 


|Vj(x,t) - u^(x,t)| = 0(h r ). 


The coefficients D , in (5.4a) can be computed directly from the 

m,k-m 

known coefficients in (5.2b) (note that B Q , k = ® k ) by successive 

differentiation of the partial differential equation and substitution — see [8] 
for detials. 

Finally, using an appropriate numerical quadrature to approximate the 
integral in (1.5a), we define ^j+ 1/2 t0 

(5.5) f j+l/2 = e k f ^ v j^ x j+l/2 , ‘ Y k T ^ ,V j+l^ X j+l/2’ Y k T ^* 

Here g^ and y^ are the coefficients of the numerical quadrature. In 
the second order case we use the mid-point rule: K = 1 , g^ = 1 , y^ = y . In 

the third and fourth order case we use the Gaussian quadrature: K = 2, 
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Bl = g >2 = Y * Yi = y (1 " 1//3), y 2 = (1 + 1//3). Note that the second- 

order accurate scheme (1.11) is identical to (5.5) with r = 2. 


Next we describe the derivation of the corrective flux g , in 

3+1/2 

(5.1b). Let R(x;w) be another reconstruction of w(x) which is at least 


r-th order accurate. Using the algorithm (3.9) we define R(x;w), its 
modified version with subcell resolution. As in the second order case (4.7a) 

A 

we define gj+j/2 * n t * ie constant coefficient case to be 


g j+l/2 = 7 l [R<x j+l/2 ' at;vn) " L(% j+l/2 " ati ' ,n)Idt 

0.6) x 

x j+l/2 „ 

= — / [RCyjv 11 ) - L(y;v n )]dy. 

X j+l/2-ar 

We note that relations (4.7b) - (4.7c) hold for any r; therefore, we can 
state that the scheme (5,1) in the constant coefficient case is identical to 
(4.2) in general. Consequently, if the reconstruction R(x;w) is exact for 

smooth polynomial data of degree r, then the ENO scheme with subcell resolu- 
tion (5.1) is exact for all initial data of discontinuous piecewise-polynomial 
functions of degree less or equal r. 

Let us denote 


?2 

(5.7a) d (y,,y 9 ) =/ [R (x;v n ) - L (x;v n )]dx, 

m l J m m 

y i 

y 2 

(5 * 7b) C m+l/2 (y l’ y 2 ) = / lR mfl (x;vn) " V XJV “ )ldx - 


Using these notations we can evaluate 


g j+l/2 


by the following expressions: 
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If 

(5. 8a) + 
unless 
(5. 8b) + 
in which 

(5. 8c) + 

If 

(5.8a)“ 

unless 

(5.8b) - 


a > 0, then 


Tg j+l/2 V X j+l/2 


aT » X j+l/2 ) » 


°j > °j-l» °j -°j+l’ F j^ X j+l/2^’ F j^ x j-l/2^ - °» 

case 

( d j (x j+l/2~ aT,x j+l/2 ) + C j+l/2 (x j+l/2 _ax ,x j+l/2 ) 

Tg j+1/2 - j if F j (x j+1/2 ) . F j ( x. +1/2 -aT) > 0 

' d j (x j+l/2" aT ’ X j+l/2 )+C j-l/2 (x j-l/2’ x j+l/2" ax) oth erwise. 
a < 0, then 


Xg j-l/2 _d j (x j-l/2’ X j-l/2 


-ax). 


j-1’ 


- °j+l* 




F j (x j-l/2 ) < °» 


in which case 


-31- 


(5.8c)~ xgj_ 1/2 


^j^ x j-l/2 ,x j-l/2 -ax ^ + ^+1 /?^ x -s-i /?“ aT » x i+i /?) 


j+1/2 j-1/2 j+1/2' 

if F j (x j+l/2 )-F j (x j-l/2' at) > 0 


' d j (x j-l/2’ x j-l/2' a,) + C 1 -l/2 <x ]-l/2’ x j-l/2" aT) other ” lse 


We observe that up to this point we have not specified L(x;w) and 
R(x;w). One possibility is to generalize the set up of the second order 
accurate scheme in section 3 as follows: Let r be the desired order of 

accuracy of the scheme. We start with a reconstruction via primitive func- 
tion R(x;w) which is one order higher, i.e., 


(5.9a) 


R (x;w) = n^Cxjw); 


here W is the primitive function of w, and H r+ ^ is the ENO interpolation 

of section 2. As before we denote the polynomial of degree r defining 

R(x;w) in I. by R.(x;w), and rewrite (5.9a) as a finite Taylor 

** J 

expansion: 


(5.9b) 


_ r D k k d k+1 

V x ;» ) ■ I e<--V' Vti Vdv 10 ' 

J k=0 J dx J 


Using (5.9b) we now define L(x;w) to be 


r-1 D, 


(5.10a) 


L (x;w) = (D + a h r D ) + £ Jf- (x - x.) 

J u r r k=Q k. 3 


(5.10b) 


a k = 


for k odd 


2 /(k+1)! for k even 


L 
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in 



(x;w) is a polynomial of degree (r - 1) which reconstructs 
to 0(h r ). We observe that 


w(x) 


(5.10c) 


,r-l 


dx 


— r L . (x;w) = D , 
r-1 j r-1 


— ■ _y w(x. ) + 0(h 2 ); 
dx r J 


this, as (4.4) does in the second order case, eliminates some of the non- 
smoothness in the reconstruction error which is due to the adaptive choice of 
stencils. Consequently, the ENO scheme based on this reconstruction is r-th 
order accurate in all Lp norms, including the maximum norm. Note that the 
correction to the first term in the RHS of (5.10a) takes care of the conserva- 
tion property (5.2a), i.e., 


(5. lOd) 


/ Lj(x;w)dx = w 


r 


Remark: There are other reasonable choices of L(x;w) and R(x;w), 


We may choose 


(5.11) 


L(x;w) = H r (x;W), R(x;w) = ^ H^xjW) 


or even 

(5.12) L(x;w) = R(x;w) ■ -g^-H^xjw); 

* 

note that the expression for g^ +1/2 
since d^(y 1 ,y 2 ) =0 in (5.7a). 


in the latter case is much simpler 
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6. EULER EQUATIONS OF GAS DYNAMICS 

In this section, we describe how to apply the scheme (5.1) to the Euler 
equations of gas dynamics for a poly tropic gas: 


( 6 . 1 a) u t + f(u) x = 0 

(6.1b) u = (p,m,E)^ 

(6.1c) f(u) = qu + (0,P,qP)^ 

( 6 . Id) P ■ (y - 1 ) (E -jpq 2 ). 

Here p,q,P and E are the density, velocity, pressure and total energy, 
respectively; m = pq is the momentum and y is the ratio of specific 

heats. 

The eigenvalues of the Jacobian matrix A(u) ■ 3f/3u are 


( 6 . 2 a) 


i^(u) = q - c, a 2 (u) = q, 83 ( 0 ) = q + 


1 /2 

where c = (yP/p) is the sound speed. 

The corresponding right-eigenvectors are 


( 6 . 2 b) 


r.(u) - ( q - c ), r 2 (u) 
\H - qc/ 


q » r~(u) 

1 2 / 

T 5 ' 


1 \ 

q + c I; 

+ q / 


here 
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(6.2c) H = (E + P)/p = C 2 /(y - 1) + j q 2 
is the enthalpy. 

To compute a left-eigenvector system {Jt^Cu)} which is bi-orthonormal 
to {r^(u)} in (6.2b), we first form the matrix T(u), the columns of which 
are the right-eigenvectors in (6.2b) 


T(u) = (r 1 (u),r 2 (u),r 3 (u)) 


and then define l^Xu) to be the k-th row in lT*(u), the inverse of T(u). 
We get 


(6. 2d) 


Jtj(u) = j (b 2 + q/ c,— b^q - l/c.bj) 

& 2 (u) = (1 - ^.bjq.-bj) 

* 3 (u) = j (b 2 - q/cj-bjq + l/c.bj) 


where 

(6.2e) bj = (y “ l)/c 2 


(6.2f) 


b« = 


1 

7 


q 2 b 


r 


Given {v^}, approximation to {u(xj,t n )}, we use (6. 2d) - (6.2f) to 

n 

evaluate the locally defined characteristic variables w^(v^) 


(6.3a) 


”i = £ k (v j )v i for 1 = J " r » 


j + r and k = 1,2,3 



- 35 - 


Note that j is fixed in (6.3a) while i varies over the points which 
are relevant to the selection of the stencil for the cell 1^. Thus, the 
eigenvector system ^k^ V j^k=l should be regarded in this context as a 

constant system of coordinates. Next we apply our scalar algorithm to each of 
the locally defined characteristic variables in (6.3a), i.e., we select a 
(possibly different) stencil for each of the characteristic variables and de- 

k 

fine Rj(x;w ) by (5.9); then we combine these scalar reconstructions by 


(6.3b) 


R (x;v n ) = l R (x;w k (v”))r, (v"). 
J ^=1 J J J 


As in section 5 we rewrite the r-th degree polynomial (6.3b) as a finite 
Taylor expansion, except that now the coefficients {D^} in (5.9) are 
vectors. With this convention in mind we proceed to define the vector recon- 
struction Lj(x;v n ) and the numerical flux by (5.10) and (5.5), 
respectively. 

We turn now to describe the vector g j+i/2* which introduces the 
subcell resolution to the numerical flux (5.1). As we have mentioned earlier 
in this paper, we use subcell resolution only in the linearly degenerate field 
(k = 2 in (6.2)) in order to improve the resolution of contact discontinui- 
ties. We do so by applying the algorithm (5.8) scalarly to the linearly 
degenerate characteristic field k = 2 as follows: 

We define 


(6.4) 


(6.5) 


F.(x) 


= R.(x.;v n )| for some k, 1 <k < r - 1, 

J dx K J J 

Z ^ 

- ^ • Z 2 (v?) {/ R j_ 1 (x;v I1 )dx + J Rj +1 (x;v n )dx - hv*} 

J X j —1/2 
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and similarly 


(6.6a) dj (y 1 ,y 2 ) - J (v” ) [R^ (x; v n ) - L^(x;v n )]dx 

y l 

y 2 

(6.6b) C nrH/2^ y l ,y 2^ = / *2^ v j ^®nH-l ^ x; ^ “ R m (x;v n )]dx for m * 


The characteristic speed of the linearly degenerate field is the flow 
velocity q; which can be of different sign in different regions. The defini- 


tion of g j+1/2 


in (5.8) is formulated as the contribution of the cell 


Ij to the numerical flux. Therefore, it is convenient to program the 
calculation of the numerical flux in two stages: First we evaluate 


(6.7) 


f = fEN 0 
r j+l/2 i+l/2 


by (5.5) for all j. Then we sweep over the mesh again and collect the con- 
tribution of each cell to the numerical flux. Using FORTRAN conventions this 
can be described by: 


If > 0 then 


(6 - 8)+ V/2- W + 7 


d j (x j+l/2 " Tq j’ x j+l/2 )r 2 (v j ); 


qj<0 


-j-1/2 


f j-1/2 “ 7 d j (x j-l/2’ X j-l/2 " Tq j )r 2 (v j ) ' 


( 6 . 8 )“ 
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Next we check whether the discontinuity condition 


(6.9) 


j > a J-l* 




J+l* 




F j (x j+l/2 ) 


is satisfied. If one or more of the inequalities in (6.9) is not true, we 
move on to the next cell. If all the inequalities in (6.9) are true we pro- 
ceed to calculate as follows 


If q” > 0, then 


(6. 10a) + 5j = 


and 


: j+l/2 (x j+l/2 " xq j» x j+l/2 ) if F j (x j+l/2 ) * F j (x j+l/2“ Tq j > 0 
c j-l/2 (x j-l/2’ x j+l/2" Tq j ) otherwise 


(6. 10b) + 


f j+1/2 f j+1/2 + 6 j r 2 (v j ) * 


If q^ < 0, then 


- 


(6.10a)" 6 -I 


and 


( X i — 1 /9 * X <4-1 ^ ^4 ^ X 44-1 /0 ^*4 ( X 4_1 /•>“ T 9 j) ^ 0 


'j+l/2 v j-1/2 1 j j+1/2 7 “ x j VA j+l/2'* r j VA j-l/2 tq j 

: j-l/2 (x j-l/2’ x j-l/2" Tq j ) otherwis e 


(6. 10b)" 


T j-l/2 T j-l/2 + 6 j r 2 (v j ) * 


Once we have completed the calculations in (6.10) we move on to the next 


cell 
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7. NUMERICAL EXPERIMENTS 

In this section, we present results of several computer experiments with 
the ENO schemes (5.5) and their modified version with subcell resolution 
(5.1); we refer to the latter as ENO/SR. 

In all these experiments we have used 


(7.1) 


a . 

3 



R(Xj ;v n ) | 


and similarly k = 1 in (6.4) for systems. In all the calculations reported 
in this section we have used a CFL number of 0.8. The continuous line in 
Figures 2-8 represents the exact solution. The circles in all Figures 
represent values of R(xj;v n ) at the time specified. 

We start with the scalar constant coefficient problem 


(7.2) ^ + U x = 0, u(x,0) = Uq(x), -1 x _< 1 

with periodic boundary conditions at x = ±1. In this case, we take 


(7.3) 


■p 

f (ui,U 2 ) = uj. 


First, 

we present numerical experiments with the highly discontinuous initial 

data 


r , ,3 2 n -1 < x < 

i -xsin(y irx ) 

1 |sin(2itx)| |x| < -j 

^2x-l-sin(3tfx)/6 < x < 1 

1 

7 

(7.4a) 

Uq (x + 0.5) = | 

• 


Note that the RHS of (7.4a) is shifted by (-0.5) for purposes of display. We 
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initialized the calculation by taking 

x j +1/2 

(7.4b) v j = h ^ _i/ 2 u o (x)dx = h^O^j+l^ “ U 0 (x j-1/2 )J 

where Uq is the primitve function of uq(x). In Figures 2 to 6, we show 
results with h = 1/30 (i.e., 60 cells) at: (a) t = 2 (after 1 period = 75 

time steps), (b) t = 8 (after 4 periods = 300 time-steps). In Figure 2, we 
show for sake of comparison the results of the MUSCL-scheme (1.11) with a 
slope s” defined by 

(7.5a) S" = m(2(Vj +1 - vp , j (v" + j- v"_j), 2(v^ - Vj_j))/h; 

here m(x,y,z) is the minmod function 

( s • min( |x|,|y|,|z|) if sgn(x) = sgn(y) = sgn(z) = s 
(7.5b) m(x,y,z) = j 

( 0 otherwise 

In Figure 3, we present results of the second-order accurate ENO scheme (4.3), 
and in Figure 4 we show results of the corresponding second order accurate 
ENO/SR (4.5) with (4.13) + . In Figure 5, we present the results of the fourth- 
order accurate ENO scheme (5.5) with r = 4, and in Figure 6 we show the 
corresponding results of the fourth-order accurate ENO/SR (5.1). 

Next we demonstrate the kind of accuracy to be expected from these 
methods in smooth problems by calculating a refinement sequence for the 
periodic constant coefficient problem (7.2) with initial data 


Uq(x) = sin(Trx) 
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In Table 1, we show the results at t = 2 with h = , -g- , -yg- (i.e., 8, 

16, and 32 cells, respectively). The quantity r c in this table is the 

"computational order of accuracy" which is evaluated from two successive cal- 

r c 

culations by assuming the error to be a constant time h ; clearly this 
definition is meaningful only for h sufficiently small. 

We turn now to present numerical experiments with the Euler equation of 
gas dynamics (6.1). In these calculations we take y = 1.4 and f R (ui,u 2 ) 
= f R0E (u 1 ,u 2 ), where 

3 

(7.6a) f R0E (u 1 ,u 2 ) tf(u 1 ) + f(u 2 ) - l 6 k |a k (u) |r k (u) ] , 

k=l 

with 

(7.6b) 6 k = £ k (u)(u 2 ~ Uj); 

here a k , £ k and r k are the eigenvalues and the left- and right- 

eigenvectors, respectively. u is a particular average of u^ and u 2 
which is defined by: 

(7.6c) q = <q/p>/</p>, H = <H/p>/</p>, c = (y - l)^^/H-l/2q Z ; 

here < > denotes arithmetic average, i.e. 

<b> = y (b i + b 2 >« 


In Figures 7 and 8, we show results of the Rieraann initial value problem 
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( U L 

x < 0 

(7.7a) 

u t + f(u) x = 0, u(x,0) = 

I, 

x > 0 


with 


(7.7b) (p L , q L , P L ) = (0.445, 0.698, 3.528); (p R ,q R ,P R ) - (0.5,0,0.571). 


These calculations were performed with 100 cells, h = 0.1, CFL = 0.8 and 85 
time-steps. In Figure 7, we show the density computed by the second order EN0 
scheme and in Figure 8 we show that of the corresponding ENO/SR. 

Finally, we present numerical solutions to the problem of two interacting 
blast waves: 


0 < x < 0. 1 
0.1 < x < 0.9 
0.9 < x < 1 

where 

(7.8b) P L - P„ - P R - 1. - q R - 0. P L - 10 3 , P M - 10- 2 . P R - 10 2 ; 

the boundaries at x = 0 and x = 1 are solid walls. This problem was 
suggested by Woodward and Colella as a test problem; we refer the reader to 
[12] where a comprehensive comparison of the performance of various schemes 
for this problem is presented. We refer the reader to [8] for a detailed 
description of the implementation of the solid wall boundary condition in the 


(7.8a) 


u(x,0) = 


ENO schemes 
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In Figure 9, we show the density at t = 0.038 calculated by the second 
order accurate ENO/SR with 800 cells and CFL - 0.8. The circles in this 
figure represent values of R(Xj,p n ); the continuous line is just the 
piecewise-linear interpolation of these values. Comparing these results to 
the solution presented by Woodward and Colella in [12], we find that it shows 
all the important features of the various interactions and thus can be con- 
sidered a "converged" solution. We use this piecewise-linear interpolation of 
the calculation with 800 cells as the "exact solution" in Figures 10 and 11. 
The circles in Figures 10 and 11 are reconstructed values of density in a 
calculation with 200 cells. In Figure 10, we show the calculation by the 
second order accurate ENO scheme; in Figue 11 we show the results of the 
corresponding ENO/SR. 

In the following, we make several remarks and observations concerning the 
numerical results presented in this section. 

(1) In all our calculations we find that the subcell resolution technique is 
capable of producing perfectly resolved linear discontinuities. Observe 
that if R(x;v n ) has a single intermediate value at a discontinuity then 
this discontinuity is perfectly resolved by R(x;v n ). 

(2) When we study the effect of higher formal order of accuracy in the calcu- 
lation of discontinuous data by the ENO schemes, we find that the most 
noticable improvement is due to the reduction in smearing of the linear 
discontinuities. However, when we compare the second order and the 
fourth order ENO/SR schemes we see that the improvement is primarily due 
to higher accuracy in the smooth part of the solution. Consequently, 
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there is not sense in going to higher order when solving a Riemann IVP. 
To justify the increased computational cost associated with higher order, 
one needs a lot of structure in the smooth part of the solution. 

(3) Comparing the solution of the interacting blast waves (7.8) by the second 
order ENO/SR to that of the PPM in [12], we find that the ENO/SR is more 
accurate. The ENO/SR highly resolves all three contact discontinuities 
in the problem, while for some reason the PPM resolves well two of the 
contact discontinuities but smears the one which results from the shock 
interaction. Another possible explanation for the difference in accuracy 
may be due to the fact that the ENO/SR is uniformly second order 
accurate, while the PPM (because of its monotonicity constraints) reduces 
to first order accuracy at points of local extremum. 

(4) The numerical results for the Euler equations of gas dynamics clearly 
demonstrate that shocks are highly resolved by the original ENO schemes, 
and subcell resolution is not needed there. In any case, the expressions 

A 

for Sj+ 1/2 to t> e modified before applying subcell resolution to a 

genuinely nonlinear field, as follows: (i) a”? should be replaced by 

the speed of the shock, (ii) Dissipation proportional to 
(aj + j/2 “ a j _ i / 2 ^ should be added to a centered rarefaction wave. 
Fortunately, if the discontinuity condition (3.9a) is met in the cell 

A ^ 

I j , then R(x;v ) is continuous at x j±j/2 an< * no interaction terms 
need be added to the numerical flux. However, one has to account for the 
fact that the wave from the interior of Ij crosses its boundaries 


during the time-step 
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APPENDIX: Derivation of the discontinuity condition 

The reconstruction R(x;w) is by definition discontinuous at 
{ x 1 / 2 } • In regions of smoothness of w(x) the jump of the reconstruction 
at x. . i /o is of the order 0(h r ). When a discontinuity of w(x) is located 


at Xj + i/2 is ot the order utn wnen a discontinuity or w^x; is located 
in the interior of I j , then the discontinuities of the reconstruction at 
X j±l/2 are fra g™ ents that w ( x )* (See Figure lb). 

In order to recover a possible discontinuity in the interior of each 
cell, we would like to associate the reconstruction with the boundaries of the 
cells { x j + 1 / 2 ^ » rat h er than the cells themselves. Let R j+i/2^ x ’ w ^ he t ' ie 
polynomial description of such a reconstruction which is valid in the neigh- 
borhood of x j+i/2* Once this is done we consider reconstructing w(x) in 

h by 


(A. 1 ) 

. _ Vl/2 <X ' , “ ) 

Rj(x;w) = • 

R j+1/2 (x; w) 

for 

X j— 1/2 1 

x < e 3 


for 

8. < X < 

X j+l/2 

if possible. 

i.e., if there is a 0^ 

such 

that 



(A. 2a) 


°* x j-l/2 — 9 j X j+l/2 


where 


(A. 2b) F(z)=I{/ 

^ Xf 


z _ x j+l/2 

/ R. L , 2 (x;w)dx + / 

X j-l/2 2 


*j+l/2 (x;w)dx} “ w j' 


If there is no such 6^, we define 
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(A.3) 


Rj(x;w) = R^(x;w). 


The only thing which is left open at this point is the definition of 
Rj + 1 / 2 (x;w). It is most natural within the framework of the ENO reconstruc- 
tion to select the "smoother" of R^(x;w) and Rj +1 (x;w), i.e.. 


(A. 4) 


Vl/2 (X;W) 


Rj(x;w) if cfj < Oj +1 


R. +1 (x;w) if >1 °j +1 


Here ck is a monotone increasing function of the "non-smoothness", such as 
(3.8a). Note that since R^(x;w) is associated with a stencil of points, 
(A. 4) is equivalent to assigning a stencil to x j+i/2> the boundary of the 
cell. 


We observe that there is a certain ambiguity in the definition of 
Rj(x;w) in (A.l) since need not be unique. We remove most of this 

ambiguity by adopting a policy of "no unnecessary changes", and agree that 


(A. 5) 


1^. i/o R * 

J-l/2 J 


or 




j+1/2 j 


=> R 


j j 


i.e., if one of the candidates for extension into Ij is no smoother than the 

A 

original R j , we just retain the original definition in R ^ . From the 

definition (A.4) of H we can rewrite (A. 5) in terms of as 


(A.5)' 


“jl'H 


or 


< Vi 


R. 

3 


R., 

3 


Hence we conclude that we define R.(x;w) 


to be (A.l) only if 
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(A.6) a.. > and a on + 1 

which is the discontinuity condition (3.8b). In this case 
and ^j+ 1/2 = R j+i and the definition ( 3 * 9 ) follows. 


Vl/2 = R j"l 



Table 1. Refinement sequence 



5.380 xlO" J 1.638 xlO J 8.800 xlO -4 5.300 xlO b 7.630 xlO 











E 



Figure la. w(x) in (3.1); circl 






















Figure 10. Second-order 
ENO. Density at t - 0.0 
with 200 cells. 
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t = 0.038 with 200 cells 
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